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The high-temperature cubic-tetragonal phase transition of pure stoichio- 
metric zirconia is studied by molecular dynamics (MD) simulations and within 
the framework of the Landau theory of phase transformations. The inter- 
atomic forces are calculated using an empirical, self-consistent, orthogonal 
tight-binding (SC-TB) model, which includes atomic polarizabilities up to 
the quadrupolar level. A first set of standard MD calculations shows that, 
on increasing temperature, one particular vibrational frequency softens. The 
temperature evolution of the free energy surfaces around the phase transi- 
tion is then studied with a second set of calculations. These combine the 
thermodynamic integration technique with constrained MD simulations. The 
results seem to support the thesis of a second-order phase transition but with 
unusual, very anharmonic behaviour above the transition temperature. 



I. INTRODUCTION 

A large class of advanced ceramics are solid solutions of zirconia (Zr02) with cubic 
stabilising oxides like Y2O3, MgO or CeO, and are generally called stabilized zirconias. 
The long list of functional applications includes high-temperature devices, thermal barriers, 
and oxygen sensors. Moreover, partially stabilized zirconias represent a new generation of 
structural materials, by far the toughest ceramic oxides, strengthened by the mechanism 
called transformation toughening. The processing and service conditions of these materials 
involve phase transformations whose underlying physics is still a subject of controversy. One 
of these is the high-temperature cubic-tetragonal phase transition which is the subject of 
the present paper. 

Zirconia is monoclinic (m) at low temperatures,liHi tetragonal (t) between 1400 and 2570 
K,0i and cubic (c) up to the melting point of 2980 K.i0 High-temperature X-ray experiments 
on stabilized zirconia revealed the existence of a c t phase transition between 2300 and 
2600 kM^ depending on the atmosphere, but the mechanism of the transformation still has 
not been fully explained. The c and t unit cells are shown in Figure |I]: note the characteristic 
tetragonal distortion of the oxygen sublattice in the t phase. 
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It is not possible to quench to low temperature the c and t forms of pure zirconia, hence 
the experiments are difficult because of the high temperatures involved. Alternatively, the 
c and t structures may be stabilized at low temperatures by impurities. The available mea- 
surements are mostly done on stabilized samples. This simplifies the experimental procedure 
but complicates the interpretation of the results, because, besides the equilibrium t phase, 
other metastable tetragonal structures are observed in stabilized crystals, denoted by t' and 
t". The former is the microstructure of a solid solution quenched from the field of stability 
of the c phase into the biphasic c + t one.S'El The t and t' forms are the same phase, they 
belong to the space group P4:2/nmc, but have different composition;0 t' is also called non 
transformable because it does not spontaneously transform to the m phase. The t" structure 
is observed in the Zr02-ErOi.5lii and Zr02-Y203lil systems, and has a cubic unit cell with 
the oxygen sublattice tetragonally distorted. 

The microstructure of samples rapidly cooled from the c-phase region presents twinned 
domains separated by antiphase boundaries. The nature and composition of these domains 
is related to the phase transition mechanism and has been a subject of controversy. Orig- 
inally they have been interpreted as the result of a diffusionless martensitic reaction.li3ll2l 
Later, Heuer and Riihle0 suggested that the transformation could be raora-martensitic: ho- 
mogeneous, massive, and displacive. Similarly, the observations of Lantieri et a/.0 were 
interpreted to mean that the c —>■ t' transformation is diffusionless but non martensitic, 
and that the transformation always goes to completion. The same authors later proposed 
that the transition could be heterogeneous of the first-order with nucleation.ll According to 
Sakuma,il the transition is instead second-order. 

The temperature evolution of the tetragonality c/a and of the anion sublattice distortion 
has been followed by Yashima et a/.Q in the Zr02-ErOi.5 system. They showed that both 
order parameters depend continuously on the temperature and suggested that "the transition 
has the nature of a higher-order phase transition" . 

Several attempts have been made in order to include this transformation in phenomeno- 
logical theories. Hillert and Sakumail expanded the free energy in terms of the defect 
concentration and assumed the transition to be second-order. Fan and Chenil used the 
time-dependent Ginzburg-Landau theory to expand the free energy of the transformation, 
treating it as a first-order one. The transformation was instead assumed to be second-order 
in the Landau energy expansion of Katamura and Sakuma.^! 

The theoretical treatment of the c ^ t transition is simpler in stoichiometric zirconia: 
this is a case similar to the c ^ t phase transition in BaTiOa, where, according to symmetry 
considerations, the transformation could be either first or second order. A free energy 
Landau expansion for zirconia, involving the tetragonality of the cell only, without the 
distortion of the oxygen sublattice, inevitably predicts a first-order transformation.!! But 
the inclusion of the latter in the Landau expansion opens the possibility for a second-order 
transition.^ As already pointed outS the coupling between the order parameters may change 
the order of the transition from second to first. 

In the case of BaTiOs it has been possible to measure the order parameters very close 
to the transition temperature0'0 and to establish the order and the mechanism of the 
transformation. Analogous experiments are difficult in pure zirconia because of the high 
transformation temperature. The neutron diffraction analysis of Aldebert and Traversei 
provides the most complete thermomechanical description of pure c and t zirconia at high 
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temperature. Aldebert and Traverse observe that: (i) The tetragonal distortion of the 
oxygen sublattice persists in the whole field of stability of the t structure, (ii) The tetragonal 
distortion of the oxygen ions vanishes in the c structure, (iii) The volume thermal expansion 
is linear and very close to isotropic up to near the transition point. As a consequence the 
c/a ratio is almost temperature-independent over a wide range of temperatures, and sharply 
decreases near the transition temperature, (iv) The isotropic Debye- Waller factors of both 
species strongly increase before the transition temperature: the authors interpret it as a 
possible structural phenomenon anticipating the phase transition which could increase the 
ionic mobility. 

The plan of the present paper is as follows. In Sec. |I| we introduce the main theoretical 
tools we have used to study the phase transition: the Landau theory of phase transforma- 
tions, the thermodynamic integration technique, the constrained dynamics and the analysis 
of the order parameter fiuctuations. The results are discussed in Sec. pT[ The phase tran- 
sition mechanism was investigated using two sets of calculations. The first one, described 
in Sec. [Ill A| , is a traditional Molecular Dynamics (MD) analysis with which we observed 
the softening of a particular vibrational frequency. The second one, described in Sec. |111 B 



combines the thermodynamic integration technique with constrained MD simulation to cal- 
culate the free energy surfaces around the phase transformation. We summarize the results 
in the final section. 



II. THEORY 
A. Landau theory of the phase transition 

1. Order parameters 

The Landau theory of phase transformationsS describes the relationship between two 
crystal structures which share a common symmetry group Go- The disappearance of a 
particular symmetry operation is quantitatively described by order parameters, which are 
zero in the high symmetry phase and become non-zero in the low symmetry one. 

Our preliminary analysisil of the c ^ t phase transformation, based on K calculations, 
showed that the transformation is driven by the distortion of the anion sublattice which 
is described by the primary order parameter 6. This is a measure of the distance between 
each oxygen atom and the corresponding centrosymmetric position it occupied in the c 
structure. The T = K calculations of certain phonon frequencies of the c phase show that 
a frequency of vibration at the X point of the BZ is imaginary.ililil This phonon, labelled 
X2 , involves the oxygen sublattice only, and is shown in Figure ||. It transforms according 
to the A2 irreducible representation of the little co-group of the X-point -D4/1. The star 
of contains three equivalent points, consequently the order parameter describing the 
tetragonal distortion has three components: S^;, Sy, and Sz- 

In transforming to the t phase, the primitive unit cell doubles, so that the phonon 
corresponding to X2 is at the F point, and is generally labelled Aig. Nevertheless, in order 
to unify the description for both the c and the t structures, we will not use this convention 
and we will always refer to the soft mode as the one, also in the t phase. 
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Besides the tetragonal distortion of the oxygen sublattice, it is necessary to capture the 
change of the unit cell shape. This is done by introducing auxiliary order parameters rji 
defined in terms of the strain tensor e. We decompose the six independent components of 
the strain tensor into irreducible representation of the Oh cubic point group in Table |I[ 

It was shown previously that, at equilibrium, each auxiliary order parameter is second- 
order in 6. From now on, the order of the expansion terms will be expressed with respect 
to the order in 6, therefore, as an example, a term like is fourth order. 



2. Energy expansion 

The Landau theory assumes that the appropriate thermodynamic potential of the crystal 
$ can be expanded in powers of the order parameters about the transition point. The Taylor 
expansion of $ must be invariant under the symmetry operations of the high symmetry 
phase. As a consequence, the allowed terms in the expansion have to be symmetry-invariants 
as well, and can be found using group theory. The terms in the energy expansion will be 
polynomials in the strains and displacements 5j of Table |. We constructed all the possible 
polynomials up to the sixth-order and symmetrized them with respect to the symmetry 
operations of the cubic point group Oh- The resulting invariants are shown in Table |T|. 

This analysis showed that all the third-order invariants are identically zero, which is a 
necessary (but not sufficient) condition for a phase transition to be second-order. We already 
mentioned that the instability appears at the boundary of the BZ, therefore it halves the 
number of symmetry elements, and this is a further condition allowing the c ^ t phase 
transition to be second-order.E2l 

In order to keep the discussion as simple as possible, at this stage we assume the phase 
transition to be second-order, truncating the Taylor expansion at the fourth-order term in 
S. The possible importance of the higher order terms will be discussed later. The energy 
expansion, expressed in terms of the basis function defined in Table ||, is as follows: 

F = F, + ^A,{6^)+j:^A,,{6^) (1) 

i=i ^ 

1=1 i=i ^ 

Fo is the energy of the high-symmetry phase and is a function of the hydrostatic strain 
rji = Tr{e). The choice of the reference volume fixes Fq and the expansion coefficients 
a2, ...,C3. In the present case, the energy F was expanded about the minimum of the 
energy-volume curve for the c structure predicted by the SC-TB calculations.il 



B. Free energy calculation 

Free energy surfaces may be calculated directly from MD simulations in terms of ensemble 
averages by using the thermodynamic integration technique.^'H Here we briefly describe how 
this method was applied to zirconia. 

The thermodynamic integration method allows us to calculate free energy differences 
between a reference state, for which the internal energy Uq is known, and another state of 
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the same system, with internal energy U. The idea is to relate the two structures with a 
switching parameter 7 which is zero in the reference state and non-zero otherwise. The free 
energy variation in the infinitesimal change d'y may be calculated using standard statistical 
mechanics: 



This is equivalent to the reversible work done for the structural modification described by 
dj, implicitly assumed to be adiabatic. By (■ ■ ■) we indicate the ensemble average, which 
has to be calculated at a constant value of 7 = 7. The free energy difference can be obtained 
by integrating the previous equation. In the general case, U{'f) is not known. The common 
strategy is to perform several constrained MD simulations at different values of 7 and then 
integrate Eq.(0) numerically. Many calculations may be necessary in order to integrate 
Eq.(^ with sufficient precision. 

A knowledge of the functional form of the energy would greatly simplify this procedure, 
reducing the number of calculations and allowing the analytic integration of Eq.(0). The 
Landau theory in combination with MD simulations can provide such useful information. 
In order to apply this formalism, it is necessary to define the thermodynamic variables of 
Eq.(|I]) from a MD run at finite temperature. Statistical mechanics allows us to calculate the 
order parameters by averaging the corresponding time- dependent ones over all the available 
atomic configurations. 

The primitive c cell is unstable with respect to 3 modes of vibration whose frequency is 
degenerate. The instability appears at the X, Y, and Z points of the primitive BZ and the 
corresponding eigenmodes distort the anion sublattice along the x, y, and z directions. In the 
following we consider a supercell which is not the primitive one, and those points, originally 
at the border of the BZ, are folded in at the F point. The eigenvectors are therefore real. 

Let us denote by u the atomic displacements from a perfect site of the high symmetry 
phase. We expand u in normal coordinates using the notation of Maradudin et alx^: 



K and label the atoms and their mass in the cell, e{K\j) is the eigenvector j at the F 
point of the BZ, and Q{j) is the corresponding normal coordinate. We will denote by \a), 
where a = x,y, z, the indexes j describing the soft modes. 

Given a general atomic configuration at time t (we now include the time in the notation), 
we define the time- dependent order parameter as the average displacement along X2 of the 
ro oxygen atoms of the cell: 



The time averages of these quantities, 6a, are the experimentally measurable order parame- 
ters which we now take as thermodynamic variables. 




(2) 




(3) 




(4) 
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The factor dU{t)/dSa entering in Eq.(^ can now be calculated at each time step by 
applying the definition of S given in Eq.(^, and by using the chain rule: 

Noting that the eigenvectors are orthonormal and that the first term of the sum is the force 
F acting on the atoms we end up with the following expression: 

i^ = v^E-F(«:,t)-e(«:|a). (6) 

Therefore the free energy gradient is calculated from the time average of the atomic forces 
projected along the Xg" mode of vibration: 



^ = v^^E-F('^,^)-eW«))_ • (7) 



Note that the above average has to be taken on an ensemble with a constant value of 
order parameter 6a, i-e. it is necessary to constrain the order parameters during the MD 
simulations. 



C. Constraining the order parameters 

The dynamics of canonical and microcanonical ensembles with fixed cell shape automat- 
ically constrain the auxiliary order parameters. On the contrary, in order to constrain the 
dynamics of the primary order parameters and then integrate Eq.(^ from the results of the 
MD simulation, it is necessary to modify the Lagrangian of the system.0 

The goal is to obtain an equation of motion describing the time evolution of a system with 
a fixed order parameter 6. This is done by extending the Lagrangian of the unconstrained 
system £": 

£ = /:«-^A,a„. (8) 

a 

The superscript u stands for unconstrained, the A's are the Lagrange multipliers to be 
calculated and the a's are the functions describing the constraints. Three of them are 
needed, one for each direction a of the tetragonal distortion: 

aait)=6ait)-5a = 0. (9) 

The Lagrangian of the constrained system is obtained from equations (^) and (|^), and the 
corresponding equations of motion are: 

MotiaiK, t) = Fa{K, t) - e{K\a). (10) 

\/ro 



Where and are the a = x,y, z components of the displacement and of the force. 
The orthonormality of the normal modes of vibration decouples the equations along the 
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three crystallographic directions, simplifying the implementation of the method. Moreover, 
since the tetragonal distortion involves the anion sublattice only, we need apply the above 
modified equation only to the tq oxygen atoms. In general, the Lagrange multipliers have 
to be found numerically, but in the present case (decoupled crystallographic directions and 
linear constraints) an analytical solution does exist. 

The expression of the Lagrange multipliers may be found as follow: (i) Advance the 
atomic positions with a fake unconstrained MD step, (ii) Use these unconstrained coordi- 
nates to find the multipliers which exactly satisfy the constraints, (iii) Use these values of 
A's to perform the true MD step which satisfies the constraining equations by construction. 
Here we specify this procedure for the leapfrog Verlet algorithm. 

Given a set of atomic positions u(t) which satisfy the constraining equations (H), the 
fake step involves solving the equation of motion corresponding to the Lagrangian By 
doing so, the set of unconstrained coordinates u"(t + At) is obtained. These are related to 
the constrained atomic positions u(t + At) as follows: 



u„(fi;,t + At) = + At) 



K{t) At^ 



e{K\a). 



(11) 



Applying the definition (Q) to these coordinates, a similar relationship may be found for the 
order parameters. 



+ At) =5l{t + At) 



K{t) Af 
roM, 



o 



(12) 



The analytic solution of the Lagrange multipliers is obtained by imposing the constraining 
equations a{t + At) = and then solving the resulting hnear equation in A: 



K{t) 



roM( 



o 



At2 



5l{t + At)-5^ 



(13) 



The substitution of (0) in (|TT]) gives the constrained coordinates at t + At in the NVE 
ensemble: 



u„(K,t + At) = <(K,t + At) 



(14) 



6l{t + At)-5^ 



e[K\a) 



It is important to notice that using this method, the expressions for the multipliers 
are functions of both the integrating scheme and any other additional constraints, such as 
thermostats. The simple case of the leapfrog Verlet algorithm described here has to be 
slightly modified in order to include the Nose-Hoover thermostat .I^^TEil The same procedure 
may be repeated for the NVT ensemble and the resulting equations of motion are: 



Ua(fi;,t) 



F«(/?,t) 



M, 



o 



(15) 



'5a,{t + At)-6a 
At2 



e K a 



where ^ is the thermostat variable. 
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D. Fluctuations 



The fluctuations of the instantaneous order parameter 6a{t) were used to calculate the 
frequency of a particular vibration directly from the MD run. The central point of this 
analysis is the calculation of the fluctuation correlation function spectrum: 

S^{u) = f e-*2-^* {6^{t = 0) 6^{t)) dt. (16) 



The above dynamic form factor is known to exhibit two important features,^ a temperature- 
dependent resonant peak at z/, and an additional central peak at z/ = 0. The relative mag- 
nitude of the two peaks depends on the transformation mechanism and on the temperature. 
This has been proved for phase transition mechanisms as different as order- disorder and 
displacive.@ Therefore, without loss of generality, following Padlewski et al^, the power 
spectrum (|16D can be modelled as a superposition of two peaks with the following functional 
form: 



CD 

-D + + [U — Ua) 

where A, B, C, D and z/^ are parameters to be fitted to the calculations. The analytical 
form of the time-dependent correlation function S{t) = {6a(t = 0) Sait)) may be found by 
substituting ( p!7D in ([16D and passing into the time domain with an inverse Fourier transform: 

Sa{t) = Ae-^' + C e-^' cos{27njJ). (18) 

The time-dependent order parameter is calculated from the MD atomic positions. The 
time correlation function of Sa{t) is then obtained using the multiple time-origin methodil 
and fitted to Eq.(|r^). The fitting procedure provides both the time correlation function and 
its Fourier transform. 



III. MD SIMULATIONS 

The polarizable self-consistent Tight-Binding mo deSi was used to perform two sets of 
MD simulations. In the first one, standard MD calculations were used to investigate the 
temperature dependence of the order parameters and to follow the softening of the mode 
of vibration up to the transition point. In the second one, we combined the constrained MD 
simulations and the thermodynamic integration technique to calculate the free energy of 
the c t phase transition, to study the nature of the order parameter fluctuations and to 
explain the high temperature stability of the c phase. 



A. Standard MD simulations 



1. Softening of a vibrational frequency 



The time evolution of a system of 96 particles with periodic boundary conditions has been 
followed at temperatures between 300 and 2200 K. The lattice parameter of the simulation 
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cell were the correspondent experimental values of Aldebert and Traverse, which, where 
necessary, have been linearly extrapolated at lower temperatures. During each MD run, the 
temperature has been constrained with a Nose- Hoover thermostat ,E3iill and the equations 
of motion have been integrated for not less then 5 ps with a typical time step of 5 fs. Near 
the transition point, the time step has been reduced to 2.5 fs and the total simulation time 
has been increased to 15 ps. 

The cell size was constrained by the relatively high number of MD runs necessary to follow 
the phase transition. In total, we simulated the time evolution of 96 particles for more than 
120 ps. As discussed later, the cell size does not change our qualitative description of the 
phase transition, and the 324-atom unit cell would have just implied a heavier computational 
effort, without adding further information to the physical picture provided by the smaller 
cell. 

We started the simulations from the crystallographic positions of the tetragonal phase 
and equilibrated the system at the temperature of 300 K. This temperature is well inside the 
field of stability of the m phase, however, during the MD simulations, the system remained 
in the t phase because of the existence of an energy barrier between the two structures. 
We calculated the vibrational frequencies of the t structure by diagonalising the dynamical 
matrix at the origin and at the borders of the Brillouin Zone (BZ) along the (100), (110), 
and (111) directions. This analysis showed that all the vibrational frequencies are real and 
that the t phase does not spontaneously distort towards the m structure. 

In this set of MD simulations we followed the approach of Padlewski et a/.El described 
in Section |II D| , focusing on the instantaneous order parameters 6a{t) [see Eq. (^)] which 
fluctuate about the mean value 6a- Figure ^ shows the typical time evolution of the primary 
order parameters for the MD run at T=700 K. Figure^ shows the fluctuation autocorrelation 
function S'(z/) and the corresponding frequency spectrum S(t) for the MD run at 700 K: the 
arrow points at z/^, the frequency which softens. It can be seen that the x and y components, 
corresponding to the transverse optical frequencies, are degenerate. 

On increasing the temperature, the softening of the frequency Uz is evident from the 
dynamic form factor, where the resonant peak shifts. At the same time, the primary or- 
der parameter decreases continuously (Figure^), as experimentally observed in the similar 
system ZrO2-12%ErOi.5.0 The calculated temperature dependence of the macroscopic or- 
der parameter 6z and of the corresponding vibrational frequency, shown in Figure |^, was 
then interpreted using the Landau theory. We found that the critical exponent for this phase 
transition is /3=0.35. According to the same theory, the critical exponent (3' for the auxiliary 
order parameters {f]2,V3) describing the tetragonality of the cell is bigger then /3. Therefore 
the c/a ratio should depend more strongly on the temperature than 6. 

As the transition temperature is approached, the decrease of the order parameter and of 
the corresponding frequency is accompanied by an increase in the order parameter fluctua- 
tions, which theoretically diverge at for a second order phase transition. As a result, it 
was not possible to follow the complete softening of the frequency: there is a temperature 
window about where, even though long MD simulations allow one to evaluate the average 
order parameters, it is not possible to calculate the frequency. In this temperature range, the 
frequency u is so low that the corresponding peak in the dynamic form factor S{i>) merges 
with the central peak and it is not possible to separate them. 

The theoretical transition temperature of f»1800 K is ~30% lower than the experimental 
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value of ^2600. □ This may be explained by noting that the first-principles calculations 
underestimate the energy difference between the c and t structures A?7*~^, which determines 
Tc.il@0 This is the ab initio energy barrier between the minima of the double well, which 
was used to parametrize the SC-TB model. In particular the SC-TB results underestimate 
the experimental AU^~'^ by 30%, which is consistent with the underestimate of the transition 
temperature. 

According to the renormalized phonon group theory,ii z/^ depends linearly on the tem- 
perature in both the regions T < and T > T^, and the correspondent slopes are related 
by the following relationship: 



However, our simulations at T > Tc suggest that the c ^ t phase transition in zirconia has 
a different behaviour from the ideal case described by Eq.(|19|), because no frequency was 
observed above Tc. 

The exploration of the high temperature region of the zirconia phase diagram has been 
carried out in two stages. As a first attempt, we continued the MD simulations on the 
system described above, simply increasing the temperature. This has been done up to 
2200 K. The time autocorrelation function (|I8]) of these simulations exponentially decayed 
without showing any structure. As a result, the central peak dominated the corresponding 
dynamic form factor, and therefore it was not possible to isolate the resonant peak at u 
from the central one. A possible explanation of this may be proposed by noticing that, 
according to Eq.(|l^), for T > Tc the slope of (^) is half that for T < T^. This means that 
the temperature window around T^, in which it is not possible to calculate the frequency, 
extends more in the high temperature field than in the low temperature one. Probably 2200 
K is still in the region of disturbance of the transition point. 

In order to verify if the frequency does eventually increase in the high temperature 
region, we studied a similar system with the same properties of that one described above 
but with a lower transition temperature. The idea is based on the following argument. It 
is well established that the relative energetics of the two phases is governed by a double 
well in the potential energy that depends on volume and c/a.&^ We studied in detail its 
dependence,^ which is also captured by the Landau expansion (|I|). Both the hydrostatic 
and tetragonal strains modify the double well in the same way: the smaller the volume (or 
the c/a ratio), the smaller the energy difference and therefore the smaller the transition 
temperature. Incidentally, this is connected to the ab initio underestimate of the energy 
barrier, which is calculated with the structural parameters corresponding to K. 

By exploiting this property of the energy surfaces, we made a new set of MD simulations 
aimed to explore the temperature range T > T^. In these calculations the volume was chosen 
to lower the transition temperature to :^1300 K and the cell was kept cubic {c/a = 1) even 
in the low temperature region. As expected, for T < T^, the linear softening of (Figure |^ 
set b) was obtained for this system as well. The slope was slightly different because in the 
previous simulations the thermal expansion of the cell was included in the description, while 
in this case the volume was fixed to the initial value. Because of this we could calculate 
the frequency up to within ^200 K of the transition temperature. The temperature was 
then increased to 2000 K. Surprisingly, even in this case, there was no structure in the 




(19) 
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autocorrelation function S{t) and the expected hardening of the frequency was not observed. 
This suggests that in the c phase the motion of the oxygen sublattice along the X2 mode 
of vibration is, in terms of S{t), uncorrelated. This behaviour will be clarified by the free 
energy surfaces described in Section [IIIB. 



B. MD simulations at constant S 



In our previous paper^3E3 we restricted the analysis of the K energy surface to one 
tetragonal invariant only. By doing so, we defined a simplified version of the energy expansion 
(0) involving the strain and one component of the primary order parameter. We then fitted 
the correspondent coefficients, 02, 041, 61, 62, Ci, and C2, to the results of total energy 
calculations. We also showed how the coupling between the primary and auxiliary order 
parameters could create a critical point where the transformation becomes first order. Here 
that analysis is extended by exploring the topology of the energy surface in the whole 5 
domain, and by following its temperature evolution through the phase transition into the field 
of stability of the c phase. This sheds light on the mechanism of the phase transformation 
and on the high-temperature stability of the c phase. 

The following results were obtained using a 12-atom unit cell with different c/a ratios (1, 
1.01, 1.02) at the K theoretical equilibrium volume of the c structure. Preliminary uncon- 
strained MD simulations were done to explore the effect of the cell size on the physical picture 
of the phase transformation described in the previous Section. Even in this small system, the 
frequency z/^ depends hnearly on the temperature and the predicted transition temperature 
is of ~ 1600 K. The effect of using a small cell is to shift to higher temperatures. This is 
consistent with the physical picture proposed in Section pi A|: the autocorrelation function 



S{t) measures the degree of correlation between the motion of the oxygen atoms along the 
X2 mode of vibration. We described how the temperature acts on S(t) by reducing the 
correlation until this is completely lost above Tc, where the corresponding frequency is soft, 
and where the structure is c. The small cell size and the periodic boundary conditions force 
the motion of atoms in adjacent cells to be correlated, and therefore counteract the effect of 
the temperature on S{t). As a result, in the small system, higher temperatures are needed 
to observe the complete softening of the frequency z/^. 

The 12-atom and 96-atom supercells have the same temperature dependence of 6z and 
i?z but the corresponding curves are shifted to different temperatures. We can therefore 
conclude that the phase transformation mechanism is the same in the two systems. We 
shall calculate the free energy of the transition from the MD simulations of the small cell, 
and assume that the resulting qualitative physical picture applies to bigger cell sizes. 

Before exploring the free energy temperature dependence, it is useful to simplify the 
complete energy expansion (|l]) by neglecting the order parameters which are unlikely to 
play an important role in the phase transition. The transformation between the c and t 
structures does not distort the cell shape as described by the order parameters {ti4,ti^,tiq). 
It is therefore reasonable to neglect them in the discussion of the following results. Moreover, 
even though the transformation between the c and the t structure does involve a change in the 
volume, the energetic contribution of the associated order parameter rji is well understood 
and has already been discussed. Apart from the K case, we will not consider the terms 
and C3 in the energy expansion. However, their possible influence on the character of the 
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phase transition in terms of softening of the corresponding elastic constant, will be discussed 
a posteriori in the final Section III B 3 . 



1. Topology of the K surface 

We start our analysis with the primary order parameter. Two sets of calculations on a 
stress-free cubic unit cell were used to fit the coefficients 02, 041, and 042. These have been 
determined by distorting the oxygen sublattice along (500) and along (550). We plot the 
resulting energy surface, which we take as starting point of our analysis, as a function of two 
tetragonal invariants in Figure In this simple case, because of the cubic cell, the three 
components of the primary order parameter are equivalent. 

The same set of calculations was then done on a tetragonal cell (c/a = 1.01), by which 
we determined the parameters 62 and C2. The latter is proportional to the elastic constant 
C . The transferability of the parameters was then checked by redoing the calculations for 
a different tetragonal cell (c/a = 1.02): Figure |] shows that the same set of coefficients fit 
the results for this cell as well. If z is the tetragonal axis, the tetragonality of the unit cell 
shortens the average interatomic distances in the transverse x, y plane and lengthens them 
along the tetragonal axis. As a consequence, the energy surface section in the transverse 
plane 5y, shown in Figure |^ (a), is similar to the reference one of Figure ^ but shallower 
and tighter, while it is deeper and broader along 5z [Figure |^ (b)]. 

Finally, the same procedure described above was used to fit the remaining coefficients hi 
and Ci by distorting the cell with respect to the order parameter r]i defined in Table |[ 

In conclusion, the static calculations show that the K energy surfaces can be captured 
by Taylor expansion up to fourth-order and are therefore completely defined by the set of 
coefficients given in Table fT|. 



2. Free energy surfaces 

The MD simulations were carried out in the temperature range from 50 K to 2000 K, 
constraining the primary and secondary order parameters. Let us first focus on the results 
for the cubic cell, commenting later on the effect of the c/a ratio. The explorations along the 
directions (500) and (550) fully determine the free energy surfaces to fourth order, therefore 
we constrained the order parameters along these directions from to 0.7 a.u., using the 
dynamics described in Section [11 C| . The quantity defined in Eq.(|^) was accumulated during 



the MD run and its time average provided the ensemble average required in Eq.(^. The 
analytical form (|l]) of the energy surface was then differentiated along the corresponding 
direction and fitted to the results of the simulations. For this particular case, the fit provided 
both the free energy gradient and the free energy itself. This is because we chose the reference 
energy as the top of the double well for a cubic crystal. The integration of Eq.(|^) provides 
the energy difference AF: 

AF = F(e,5)-F(e,5 = 0) = (^^^dS. (20) 
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We arbitrarily set to zero the integration constant for the cubic cell F{e = 0,5 = 0). If 
the cell is tetragonal the same constant is C2C2(e^)/2. Therefore, the fit of the free energy 
gradient to the MD simulations of the cubic cell provided the coefficients 02, 041, 042 and 62 
at the corresponding temperature. 

The fit to the computed results along the (500) and the corresponding free energy pro- 
files obtained from Eq.(pOD are shown in Figures § and ^ respectively. The corresponding 
expansion coefficients are included in Table Also at high temperature the fourth order 
energy expansion well describes the calculated free energy gradients. The temperature acts 
on the double well by gradually reducing the energy difference between the distorted and 
undistorted structure. As expected, the energy surface is very flat near the critical temper- 
ature, but, surprisingly, it remains quite flat even at higher temperatures, well inside the 
fleld of stability of the cubic phase. The energy surfaces above the transition temperature 
are highly anharmonic. 

The unconstrained MD simulations described in the previous Section generated the tem- 
perature dependence of the vibration frequency (Figure |^) up to the transition point 
only. Above the critical temperature it was impossible to calculate the frequency Vz-, the 
dynamic form factor S{u) being dominated by a wide central peak. The soft frequency Vz, 
together with the large fluctuations of the order parameters 5^, suggests a disordered dy- 
namics of the oxygen sublattice along the mode of vibration. The anharmonicity of the 
energy surfaces explains this behaviour. It may also explain the fact that the optical-phonon 
branches are not experimentally observedlll in cubic stabilized zirconia. 

As in other perfect fluorite structures,ii the vibrational motion of the anions in c zirconia 
appears to be anharmonic. This is consistent with the neutron powder diffraction experi- 
ments of Kisi and Yuxiangll on cubic stabilized zirconia (Zr02-9.4%Y203). They measured 
the temperature dependence of the Debye- Waller factor and proposed different models to flt 
the data. Both a simple Debye model and a Debye model plus a static disorder component 
provided a poor flt of the data. A radical improvement of the flt was obtained when an 
isotropic anharmonic vibration of both species was included in the description. 

We may further investigate the nature of the double well temperature dependence by 
splitting the free energy into its energetic and entropic contributions. The time average of 



the internal energy U from each constrained MD simulations is plotted in Figure 11. It is 



clear that the double well in the internal energy is present even at T > Tc. The internal 
energy double well is relatively insensitive to the temperature and the K coefficients of 
Table ^ provide an excellent fit of the finite temperature results for both the low and high 
temperature structures. From the definition of F, the difference between the free energy 
calculated with Eq.(PD[) and the internal energy gives the entropic contribution TS plotted 



in Figure 12. The high temperature stability of the c phase is therefore ensured by this 



entropic term which change the shape of the energy surface from double to single welled. 



3. Coupling to the elastic strains 

The flt of the Landau energy expansion (|lD to the calculation results has allowed us to 
follow the temperature evolution of the free energy surface through the phase transition. 



The fltting coefficients for each temperature are included in Table III. We can see that at 



the critical temperature of ~ 1600 K, when a2 goes to zero, the fourth order coefficient 04 
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is positive. Therefore, in this cell, the phase transition is diffusionless displacive and second 
order. 

If the thermodynamic potential F is expanded in terms of elastic strains as well as 5, the 
coupling between thenrimary order parameter and the strain in renormalizes the fourth- 
order coefficientilSit and could make it small or negative near the transition point. As a 
result, the transformation may become first-order. This could happen if the temperature 
reduces an elastic constant q. 

In order to see this, let us consider a tetragonal cell at the reference volume, whose 
tetragonal axis is z and where the oxygen sublattice is distorted along the order parameter 
direction (OO^^). With these restrictions the energy expansion (0) has a simplified form: 

F = F, + ^61 + ^6t + h 251 m+'-lvl + O (sf) . (21) 

At equilibrium the two order parameters 5z and 772 are not independent and the relationship 
between the two may be found by imposing the equilibrium condition: 

= ^ r/2 = '-5l 22 

dr]2 C2 

A similar procedure may be repeated for the hydrostatic strain exx + (-yy + Czz = ^^id the 
substitution of Eq.(^) back in Eq.(^), together with the corresponding one for r^i, clarifies 
the combined effect of the coefficients hi and q on the renormalization of the fourth order 
coefficient: 

f = fo + ?«+f^-^-^l*: + 0(a. (23) 



It may be verified that the above formulation is independent of the initial choice of the 
tetragonal axis. 

The coefficients ci and C2 are proportional to the bulk modulus and to the elastic constant 
C", respectively. If one or both of these quantities significantly decrease with temperature, 
the correspondent term in Eq.(|23|) may dominate the sign of the fourth-order coefficient, 
making it very small or negative. At K the coupling terms 2h\/c2 and h\/2ci reduce the 
fourth-order coefficient by 16% and 4% respectively. Because of this difference we focused 
our attention on the elastic constant C . Experimentally, the high temperature data§'iiil 
do not show any anomalous temperature dependence of the elastic constants, with a general 
decrease of ~ 15 — 20% between 300 and 1700 K. If this is valid for pure zirconia as well, 
we may anticipate that the degree of softening of the elastic constants will not affect the 
character of the phase transition. 

The calculations described above for the cubic cell have been repeated for two tetragonal 
cells with c/a = 1.01 and 1.02. The constrained MD simulations along the (500) direction 
allowed the fit of the coupling coefficient 62 to each temperature. Therefore it has been pos- 
sible to obtain the temperature and c/a dependent free energy curve along this direction. 
However, in this case we are interested in the relative position of the free energy surfaces for 
the different cells, depending on the integration constant F(e, 5 = 0) = C2 C2(e^)/2. In prin- 
ciple one could calculate C2 by a simple thermodynamic integration, but this would require 
monitoring the stress, which is not yet implemented in the current program. Therefore we 
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decided to continue the analysis by choosing an extreme scenario, which would be a strong 
temperature dependence of C . With a large safety margin with respect to the experimental 
values, we linearly reduced this elastic constant by 75 % from 200 to 50 MPa between and 
2000 K. 

With this assumption, the sections of the free energy surfaces at different temperatures 
between 500 and 2000 K are plotted in Figure |13] as a function of 5 and c/a. As experimen- 
tally observed,!^ the tetragonal cell with c/a = 1.02 is thermodynamically stable up to very 
near where the minimum configuration goes from (c/a = 1.02, 5 = 0.4) to (c/a = 1,6 = 0) 
quickly but continuously. We believe that the sudden change of the order parameters, due 
to the fiat energy surfaces near the transition point and therefore to the anharmonicity of 
the material, may explain the fact that this phase transition has been considered to be first 
order in the early studies. 

It is interesting to note that both the a^i and the coupling coefficient 62 decrease with 
temperature (Table |IT|). Because of this, even with the large postulated softening of C, the 
renormalization of the fourth-order coefficient in Eq. (p3D does not increase with temperature. 
From the data of Table fT|, we find that at Tc, the term 2h\/c2 is still only 10% of a4i/4, less 
than in the K case. These results show that even a large softening of the elastic constant 
C does not change the character of the phase transition, which remains displacive second 
order. 

Figure |13| shows also that, above the transition temperature, the minimum energy cor- 
responds to the cubic cell and therefore the results about the high temperature structural 
stability of the c phase discussed in the previous Section remain valid. 



IV. CONCLUSIONS 

The c ^ t phase transition of pure stoichiometric zirconia has been studied from differ- 
ent theoretical perspectives. Both symmetry arguments and the lattice dynamical analysis 
suggested that this transformation might be second-order. In the K perfect t structure, the 
primary order parameter S has a clear definition in terms of the displacement of the oxygen 
atoms away from their centrosymmetric position, which define a zone boundary phonon Xg". 
But the same definition cannot be directly applied to a finite temperature atomic config- 
uration. Instead, we defined the more general macroscopic thermodynamic variable as an 
ensemble average of the displacements of the O atoms projected onto the X2 zone boundary 
phonon coordinate. The corresponding frequency of vibration was then calculated directly 
from the MD simulations by applying standard statistical mechanics techniques. 

The temperature evolution of both the equilibrium order parameter 6 and of the corre- 
sponding frequency u was then followed during the MD run. The results of these simulations 
have been interpreted with the Landau theory: the critical exponent of 0.35 fitted to the 
results is between the extreme values of 0.5 and 0.25 corresponding to a second-order and 
to a tricritical phase transition. 

Approaching Tc from the field of stability of the t phase, the order parameter 6 grad- 
ually decreases up to very near the transition point, as in a displacive second-order phase 
transformation. In a temperature window about Tc, the large fluctuations of the order pa- 
rameter degrade the quality of the averages achievable from MD simulations, but we were 
able to observe z/^ decreasing linearly with T by some 70 %. In contrast to the prediction 
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of the Landau theory, no increase of that frequency in the field of stabihty of the c phase 
was observed. At high temperatures the dynamics of the oxygen sublattice revealed a high 
degree of mobility of the anions and a low correlation between their motion along the 
vibrational mode. 

In order to clarify these observations, we calculated the free energy surfaces relating 
the two structures at different temperatures, combining constrained MD simulations and 
the thermodynamic integration technique. These calculations showed that the high tem- 
perature stability of the c structure is due to the entropic contribution TAS and not to a 
variation of the internal energy profile. Moreover, we showed that the energy surfaces of the 
c phase are highly anharmonic, not only at the transition temperature, but also well above 
Tc. This confirms and may explain in terms of thermodynamic quantities the absence of 
the optical modes of vibration in the experimental spectra of Liu et a/.,0 and the postu- 
lated uncorrelated "fluidlike" motion of the oxygens about their centrosymmetric position. 
Similarly, the experimentally observed increase of the ionic conductivity and the "structural 
phenomenon" mentioned by Aldebert and Traverse! may possibly be connected with the 
soft dynamics of the oxygen sublattice caused by the fiat energy surfaces. 

Our analysis revealed the peculiar character of the c ^ t transformation in zirconia. 
Approaching from below, the variation of the free energy surfaces seems to support the 
thesis of a t — > c displacive second-order phase transition. On the contrary, approaching 
Tc from above, it is not possible to follow the softening of a particular phonon mode. The 
entropic term TAS" eliminates the double well in the free energy, and frees the atoms to move 
along that mode without an energy cost and without a well defined vibrational frequency. 
This is more akin to the high temperature phase in an order- disorder phase transformation. 
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TABLE L Order parameters for the c -fr^ t phase transition decomposed into irreducible repre- 
sentations of the Oh cubic group. 
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TABLE II. Polynomials in the order parameters of Table | which are invariants under the set of 
transformations belonging to O^. 
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TABLE III. Coefficients for the Landau energy expansion (|T|) at the K equilibrium volume 
of the c cell: 02 and 042 in Ry/ag, 041 in Ry/ag, and C2 in Ry, where ao is the Bohr radius. 
The coefficients hi = —0.363 Ry/ao and ci = 16.768 Ry complete the K set. See text for the 
temperature dependence of C2. 
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FIG. 1. Cubic and tetragonal structures of Zr02- Light and dark circles denote oxygen and 
zirconium atoms respectively. Arrows represent the structural instability of the oxygen sublattice 
along the mode of vibration. 




FIG. 2. Time-dependent order parameters at 700 K: 5x and 5y oscillate around 0, 6z around Sz, 
the value of the macroscopic order parameter. 
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FIG. 3. Time correlation functions (top) and corresponding Fourier transforms (bottom) of the 
time-dependent order parameters 6x, dy, and 6z- 
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FIG. 4. Temperature dependence of the macroscopic order parameter S^- The symbols (•) are 
the results of the calculations. The continuous solid eye-line is extrapolated in the region near 
where the large fluctuations in Sz make the averaging procedure inaccurate. 
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FIG. 5. Calculation results of the frequency squared p| vs. temperature for two simulation sets 
: (a) tetragonal cell with temperature-dependent lattice parameters taken from experiment,! (b) 
Cubic cell with temperature- independent lattice parameter (see text). 
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FIG. 6. Section of the K energy surface for the cubic cell. The isoenergetic contours are plotted 
on the base every 0.3 mRy/Zr02. 




FIG. 7. K energy surfaces for a tetragonal cell with the tetragonal axis along z: (a) section in 
the 6y,6x plane, (b) section in the 6z,6x plane. The isoenergetic contours are plotted on the base 
every 0.3 mRy/Zr02. 
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FIG. 8. Transferability of the K energy-expansion coefficients between different tetragonal 
cells with the tetragonal axis along z. Projections of the corresponding energy surfaces along the 
high-symmetry order-parameter directions (00(5) (a), and ((550) (b). 
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FIG. 9. Temperature-dependent free energy gradients calculated using Eq.(^) and corresponding 
fit via the analytic form derived from the Landau theory (||). Projection along the order parameter 
direction (005). 
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FIG. 10. Temperature evolution of the free energy profiles projected along the order parameter 
direction (005). The symbols (•) are the K calculations, the thick solid lines are the result of the 
thermodynamic integration and correspond to the temperatures 50, 500, 1000, 1500, and 2000 K. 
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FIG. 11. Double well in the internal energy along the order parameter direction (005): the solid 
line correspond to the K calculation, the symbols are the averaged internal energies from the 
constrained MD simulation. 
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FIG. 12. Entropic contribution to the phase transition obtained by applying the definition of 
the Helmholtz free energy AF = AU — TAS to the data shown in Figures |10| and 0. 
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FIG. 13. Free energy profiles projected along the order parameter direction (00(5) for cubic and 
tetragonal cells below and above the transition temperature. 
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